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The method of smoothed particle hydrodynamics (SPH) is developped appropriately for the 
study of relativistic heavy ion collision processes. In order to describe the flow of a high energy but 
low baryon number density fluid, the entropy is taken as the SPH base. We formulate the method 
in terms of the variational principle. Several examples show that the method is very promising for 
the study of hadronic flow in RHIC physics. 
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Hydrodynamic descriptions of high energy hadronic and nuclear collisions have a rather long history [EJ . Although, 



I. INTRODUCTION 

o : 
o 

• from theoretical point of view, it is not a trivial matter to justify their validity, they have been successful in reproducing 
certain features of these processes, such as the energy dependence of the average multiplicity and the transverse-energy 
distributions. More recently, relativistic fluid dynamics have become an important tool for the analysis of relativistic 
heavy- ion-collision processes (for example, and references there in). In these processes, the nuclear matter is 

expected to be compressed and heated up close to those states of the matter realized in the Big Bang Era of the 
Universe. Some laboratory data of these fascinating processes have already been obtained in the series of CERN 
experiments [g] , and data at even higher temperatures and densities will soon be available in the forthcoming RHIC 
■ experiments. 

^ ' The relativistic hydrodynamics is a description based on local conservation laws, together with the hypothesis of 
ON ■ local thermodynamical equilibrium. The conservation laws are written in terms of the four-divergence of the energy- 
momentum tensor. The resulting system of equations is highly nonlinear and analytical solutions are only available for 
some very particular and limited configurations and equations of state Thus numerical approaches are resorted 

to but they usually depend on some sophisticated techniques specific to some symmetry involved in the problem. 
When no symmetry is involved, these methods become computationally very expensive. However this is exactly the 
case when we are challenged to a realistic three dimensional simulation for nuclear-collision processes. There, we 
expect no geometrical symmetry so that a full 3D calculation is required. 

One basic point in the hydrodynamic approach of relativistic nuclear collisions is that its principal ingredients, i.e., 
the equation of state of the matter and the initial conditions for the dynamics are not quite well known. On the 
Oh, contrary, we apply the hydrodynamic models to infer these informations on the properties of the matter in such a 
highly condensed and excited state. Thus we need to perform many hydrodynamic-model calculations for different 
equations of state and initial conditions to compare with the experimental data. In such a process, we actually don't 
need the very precise solution of hydrodynamic equations, but a general flow pattern which characterizes the final 
configuration of the system as a response to a given set of equation of state and initial conditions. We are not interested 
in, and probably neither able to analyze at least in the present stage, any very precise local feature (for example sound 
ripples, small local perturbations, etc.) in these models. Extremely local properties, therefore, should be averaged out 
without spoiling the general flow pattern. This is particularly the case when one considers the questionable validity 
of the local thermal equilibrium in the problem of interest here. In short, for the study of hydrodynamic models of 
relativistic nuclear collisions, we prefer a rather simple scheme of solving hydrodynamic equations, not unnecessarily 
too precise but robust enough to deal with any kind of geometry. From this point of view, it is stressed in Ref. Q the 
advantage of variational approach to the relativistic hydrodynamics. 

Among many numerical approaches, the Smoothed Particle Hydrodynamics (SPH) |^,|9| fits perfectly in with the 
variational formalism. As a consequence, this approach presents many desirable profiles of the variational approach, 
such as simplicity and robustness with respect to the changes of geometry, as well as the possibility of smoothing 
out undesirable local degrees of freedom. Furthermore, the SPH parametrizes the matter flow in terms of discrete 
Lagrangian coordinates (called "particles") attached to some conserved quantity, such as baryon number. In this 
aspect, the SPH method is the best to the studies of relativistic nuclear collisions, where a extremely compressed and 
high-temperature hadronic matter expands into a very large space region. This is one of the principal advantages of 
the SPH over the other space-fixed grid algorithms when applied to the study of RHIC physics. The SPH algorithm 
was first introduced for astrophysical applications and used by now in several calculations, such as fragmentation 
of asteroids, supernova explosions, collision of neutron stars, etc. Several studies of the SPH algorithm for the 
ultrarelativistic regime of hydrodynamics have been done M-Of . 



1 



However, some specific aspects of the relativistic heavy ion collision processes deserve attention when applying 
the SPH. One of them is that, in the ultrarelativistic regime of central collisions, a large fraction of the incident 
energy is converted into produced particles. In particular, in the mid-rapidity region of central collisions, the most 
of these energies are in the form of produced pions and only a very small portion is carried by baryons. This is a 
very unfavorable situation to the conventional formulation of the SPH algorithm, where particles are associated to 
the matter defined by the conserved quantity, such as the baryon number. A direct application of the SPH based on 
the conserved baryon number may fail in the null baryon-number, pion-dominated region. 

As mentioned above, the basic point of the SPH method is to introduce a set of "SPH - particles" which follow the 
flow of the fluid. However, the definition of a flow does not necessarily rely on the conserved quantity. For example, 
according to Landau []l8f a flow is defined in terms of the local Lorentz frame in which the energy-momentum tensor 
becomes diagonal. In this paper, we explore this aspect and formulate the SPH in terms of any extensive quantities 
defined in the Landau comoving local frame. We derive the relativistic SPH equations using the variational principle 
, taking the matter flow as the variable. We argue that the quantity we attribute to the SPH particles convenient 
for relativistic heavy-ion collisions are the entropy of the fluid. In this way, we can follow directly the entropy content 
and its change due to the dissipation mechanism, for example, the shock wave. This is particularly interesting for the 
system where the first order phase transition is present |14|~|lq] . 

Another specific aspect of relativistic heavy-ion collision processes is how to set the initial conditions for the hydro- 
dynamic motion. As mentioned before, this is not a solved problem. However, generally speaking, the hydrodynamic 
regime will be established from its initial excited state of the microscopic QCD degrees of freedom, only after a certain 
relaxation time r re i ax . This means that the hypersurface of the onset of the hydrodynamic regime is characterized by 
a nearly constant local proper time rather than by a constant global coordinate time t. For ultrarelativistic regime, 
the difference becomes crucial. In the limit of very large initial energy, the Bjorken scaling solution is expected to 
be a good first order approximation at least for the longitudinal motion. Therefore, we shall write the hydrodynamic 
equations in terms of the Bjorken scaling coordinates r = y/t 2 — z 2 and r] = 1/2 \n(t + z)/(t — z) , replacing t and z, 
as far as the longitudinal motion is concerned. These variables are also suitable for the setup of the initial conditions. 
Here, we show that the variational approach is also useful for the derivation of the SPH equations for an arbitrary 
curvilinear coordinate system. 

We organize the present paper as follows. In Sec. II, we introduce the relativistic variational formulation with the 
SPH parametrization of one of the thermodynamical extensive quantities. In particular, we argue that the entropy 
is a convenient extensive quantity for the application of relativistic heavy-ion collisions. In this case, the change of 
entropy due to dissipative processes should be taken into account. This is particularly important in the presence of 
a shock wave or the first order phase transition. For this purpose, we deduce the entropy based SPH equations in 
the presence of bulk viscosity. In Sec. Ill, we apply our formulation for several relativistic systems such as one- and 
three-dimensional Landau model, the one presenting longitudinal Bjorken scaling solution plus transverse expansion, 
and one-dimensional relativistic fluid with shock waves. These results are compared with known solutions. Sec. IV is 
dedicated to the conclusions and perspectives of the present work. 



II. VARIATIONAL DERIVATION OF THE RELATIVISTIC SPH EQUATION 



A. SPH Representation of Extensive Variables 



The non-relativistic SPH can be formulated in terms of the variational principle. We show here that the relativistic 
SPH can also be formulated in this way 0. This guarantees that the SPH coordinates {r a (i)}are the optimal dy- 
namical parameters to minimize the SPH- model action. The other advantage of this approach is that it automatically 
leads to the symmetrized form of the SPH equation to conserve linear and angular momenta of the system JT^ . This 
is the natural consequence of the Lorentz scalar nature of the action. 

In the relativistic hydrodynamics, we assume that at any space-time point x = jr, t} there exists a local reference 
frame where the energy-momentum tensor becomes diagonal and takes the form |18| 



T» v (r, t) = 



( e \ 

P 

P 

Vo P) 



(1) 



where e and P are the (proper) energy density and the pressure of the fluid. From the hypothesis of local equilibrium, 
we assume that the thermodynamical relations are valid in each local frame. 
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Let A be an arbitrary thermodynamical extensive quantity of the fluid such as baryon number, entropy or specific 
volume. The amount of A contained in the infinitesimal volume element dV is denoted by dA so that the corresponding 
(proper) density a is 

dA 

(2) 

This is related to the corresponding density a* measured in the space-fixed (calculational) frame as 

a* (r, t) = ja, (3) 

where 7 is the local Lorentz factor associated to the flow. 

In the SPH representation, we parametrize this density by the following ansatz, 

a*(r,t) = J2 v i W (^-n(t);h), (4) 

i 

where W (r — f';h) is a positive definite kernel function peaked at r = r ' with the normalization 

d 3 rW {r-r'-h) = 1. (5) 



The parameter h represents the width of the kernel. In the limit, h — > 0, we have 

lim W (r- r '; h) = S 3 f ') . 

h— *0 

As will be seen later, it is convenient to choose an even function for W. It can be a Gaussian function in x — \f— r'\ /h, 
but in practice we often use the B-spline functions |17[ . The role of W with a finite value of h is to introduce a sort 
of short-wave length cut filter in the Fourier representation of the density a* . The total amount of A of the system is 
obtained by integrating Eq.(||) over the whole space, 

Atot = f d 3 fa* (r, t) =J^ v i ( 6 ) 

Physically speaking, it is clear from the above expression that we are replacing the system of continuous fluid by a 
collection of "SPH particles" each of which carries portion of the extensive quantity A. 

The velocity of these particles are identified as the velocity of the fluid at their position ?i(t), 



so that the Lorentz factor of the i-ih particle is given by ji = — vf- If A is a conserved quantity then {z/j} 

should be constant in time. Then 



da* 
dt 



{r~n{t);h) 

i 

= -V ■^v i v i W{r-r i {t);h) 

i 

= -v • j A , 

where 

i 

is the current density of A. Thus, the continuity equation is satisfied by the ansatz, (f|), together with Eq.(Q). If A is 
not conserved, then v[s are not constant in time. In this case, the continuity equation has the contribution from the 
time derivative of v[s as 

^+^-jA = J2hW(f-n(t);h). (8) 
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We consider the set of time dependent variables {fi,i— 1, ...,n} as the variational dynamical variables and their 
equations of motion are determined by minimizing the action for the hydrodynamic system. Here, {z^} are not 
dynamical variables and are determined by the initital conditions together with the constraints for the variational 
procedure (see the later discussion in 2.3.). 

It may seem that the spatial extension of the particle i is identified to that of the kernel W (r — f^; h). However, 
the density of A at the position of the particle i is given by 

a^Y^VjWW-fj), (9) 

3 

from Eq.(Q). Therefore, we may define the "specific volume" V% of the extensive quantity A associated to the particle 

i as 

v< = - = r*Z~ -v (10) 

a i Isj V 3 W \ r i - r j) 

Any other extensive quantities carried by the particle i can be calculated easily. Let o* be the density of another 
extensive quantity, say O, measured in the calculational frame. Then the amount of O carried by the particle i is 

Vi {0/A) i = v, {o/a) t = Vi (o/a)* , (11) 

so that the density distribution of O in the calculational frame is expressed as 

o* (r, t) - 5>< (o/a)* W (f - n {£)) . (12) 



B. SPH Action 



The relativistic hydrodynamic equations can be obtained by the variational principle for the action |Q 

I = - [ d 4 xe, (13) 



with respect to the matter density distribution, subjected to the constraint expressed by the continuity equation. 
Here, for the sake of simplicity, we discuss first the case of the Minkowsky metric. The argument can readily be 
generalized for more general coordinate systems (see Sec. III). The Lagrangian of the system is 

L= - J d 3 re. (14) 

We may consider e as the Lagrangian density in the space fixed frame. Therefore, the corresponding Lagrangian for 
the system of SPH particles can be taken to be 



Lsph ({^'^}) =-J2 v i( £ / a *)i 



-?(?).• 

where {r^ (t)} are the dynamical variables and Ei = vi (s/a) i is the "rest energy" of the particle i. The SPH model 
action is then 



iSPH 
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C. Variational Procedure 



The Lagrangian, given by Eq.(15), seems to show that the system is as if it is just a sum of free particles of "rest 
mass" Ei. However, there exists a basic difference. Namely, E^s are not constant with respect to the variation of 
positions. This is because the energy contained in each particle is determined by the configuration of other particles 
through thcrmodynamical relations. More explicitly, variations in the particles positions {fi, i = 1, n} cause change 
of the volume occupied by each particle, which in turn modifies its energy (rest energy). 

For our variational procedure, we consider changes of quantities associated only kinematically with virtual variations 
of the configuration. That is, except for the energy and the volume, all extensive quantities such as the entropy and 
the particle number should be kept constant while variations of {fi} are taken. This implies that for variations of fi, 
we keep f j constant if A is not the energy or the volume. Now, if we take A as the volume, then its changes associated 
to variations of fi are described by Eq.(^) where Vi should be understood as the initial volume of the particle i. 
Therefore, even for the case A is the volume, vi should be kept constant for variations of fj. In short, except for the 
case of the energy, the parameters Vi should be kept constant while variations of {fi} are taken. If we take A as the 
energy, then the change of Vi should be calculated as a function of the variation in {fi}. This introduces an additional 
complication without any practical merit. Therefore, for the sake of simplicity, we restrict ourselves to the case where 
A is not the energy but any other extensive quantity. 

We can write the change of energy associated with a virtual change of volume 8V as 



SE = -P5V + SW 



(17) 
(18) 



where SW is an additional work to change irreversibly the volume and P e f / is the effective pressure. If the change of 
volume is performed in a quasi-static way and there exists no dissipative force, then this effective pressure coincides 
with the usual pressure P. However, if there exists some irreversible process associated with the volume change, then 
P e f f ^ P and we may write 



eff 



P + Q, 



(19) 



where QSV is the energy change due to some irreversible process. For a quasi-static adiabatic process, Q = 0. The 
introduction of Q is important when we deal with shock phenomena (see Sec. III). 
The variation of volume SV for constant iVs is calculated from Eq.(|l(i|) as 




^E^ 



(20) 



where Wij = W (fi — ff, h). Using these relations, the variation of the action (16) with respect to {n,i — 1, n} 
becomes 



SISPH = - I ' *X) ~ f" ( P + Wi - E > — 



■/*!>•{- 



e + P + Q 



a 



liVi 



i (P + Q\ J_ fP + Q 

' 2 

i Ij 



2 I 9 I 9 \ 9 

7i V a Ji 7," V a 



Here, we used the symmetry of W and the property 

ViWij - V,ll, ; . 

The requirement SIsph = for Vi5r a leads to 



(21) 



d_ 

dt 



P + Q 



■E 



% 3 



P + Q 



P + Q 



(22) 
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The corresponding hydrodynamic equation in the continuum limit is 

d»T» v = (23) 

where 

= Q [u^u v - g" v ] (24) 



is the stress tensor for the bulk viscosity Q. This equation is the same as those discussed in Refs.( Jl l| — |l 3[| ) - 

Equation (p^ ) should be complemented by an equation which expresses the conservation of energy (and other 
thcrmodynamical quantities, if any diffusion process is present). The energy conservation is written as 

On the other hand, from the second law of thermodynamics, we should have 

^i = - P ^i +T .^i +fjLi ^i (26) 

dt dt dt dt 

in equilibrium, where T,S,/i,N are the temperature, entropy, chemical potential and particle number, respectively. 
Thus, we get 

dS, dNi dVi 

T ^ + ^ = - Q *nr- (27) 

In the absence of particle diffusion, the chemical equilibrium requires 

dN 

,f=0 (28) 

and, in this case, 

D. Entropy Representation of SPH Equations 

As we have mentioned in the Introduction, in applications to ultrarelativisitc nuclear collisions the baryon number 
is not a suitable quantity to represent the hydrodynamic flow, since most of the energy content is in the form of non- 
baryonic matter. This is particularly so in the central rapidity region. We may consider the energy content itself as 
the SPH base. However, as mentioned before, this choice introduces an additional constraint between the coordinates 
{fi} and the extensive parameters {vi\ of SPH particles, due to the energy conservation, and not desirable from the 
practical point of view. Therefore, we propose to take the entropy as the suitable extensive quantity for the SPH 
representation. 

Let s* (t,r) be the entropy density in the space-fixed (calculational) frame. From Eq.(^), we have 

s* = s* (t,n(t)) 

3 

and the equations of motion for {r*j} are given by Eq. (p2|) , which we write in the form, 

§ = (30) 



dt 

where 



ViWij, (31) 



~, = vi 7,^ (32) 



G 



is the momentum density associated with the i — th particle. 

For the variational procedure, i/j's are kept constant but it does not mean that they are constant in time. When 
there exists some non-adiabatic process, then Q^O and we have to use the energy conservation equation, Eq.^) to 
determine the change of these parameters {vi}- We have 



1 dvj _ Qj 
V t ~dt ~ ~T^* 



(33) 



where 



- ldVi 
1 ~ V, dt 

= d^, (34) 

Eqs. (30,31,^) constitute a system of first-order differential equations for r*;,7r.; and where the velocity Uj should 
be determined algebraically from Eq.([32|). Further, we have to specify the equation of state, for example, 

E = E(N,S), 

and the dissipation pressure, 

Q = Q(N, S,0). 
We will discuss these points later in practical examples. 

E. SPH Equation for Generalized Coordinate System 

The variational procedure can readily be extended to coordinate systems with non-Cartesian metric. The use of the 
generalized coordinate system is particularly important when we consider realistic initial conditions for simulations 
of RHIC processes. As we know, in a relativistic heavy-ion collisional process, the initial state is a cold, quantum 
nuclear matters flying separately. Just after the collision, the hadronic matter stays at a highly excited state and 
the materialization occurs only after ~l/m/c in the proper time. Therefore, the local thermodynamical state would 
emerge for some local proper time and not for the global space-fixed time t. Thus, it is important to choose a 
convenient coordinate system for the description of relativistic heavy-ion collisions. For example, it is often used the 
hyperbolic time and longitudinal coordinates as is described later. 

Let us consider a more general coordinate system, 

ds 2 = g liv dx ll dx v . (35) 

However, in order to unambiguously define the conserved quantities, we consider only the case when the time-like 
coordinate is orthogonal to the space-like coordinates, 

S„o = 0. (36) 
The action principle for the relativistic fluid motion can be written as 

j4 



SI = -5 J d 4 x^g~£ = , (37) 
together with the constraint for the conserved entropy current, 

= {V=g^) = o (38) 
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±=d T {^—g S1 ) + -L ^2 d i (V^gsjv 1 ) = , (39) 
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where 



(40) 



and we use the notation, 



t = x , 7 = u . 



The generalized gamma factor 7 is related to the velocity v a through u^u^ = 1, so that 

1 



7 = 



V500 - ^gv 



(41) 



where — g is the 3x3 space part of the metric tensor. That is 

(^) = ("S° -°g)- ( 42 ) 

Let us now introduce the SPH representation. We may, for example, express the entropy density by the ansatz 

V^sj = s* -> s* SPH = ViW (r - n (t)) , (43) 

i 

or by 

s 7 = s* 4 Pff = ^ i/jW (r* - r* (r)) , (44) 

as well. These two possibilities, besides others, are simply different ways to parametrize a variational ansatz in terms 
of a linear combination of given functions W (r — f j (r) ) . The most important property of an ansatz should be that 
W satisfies the normalization condition imposed by the basic conserved quantity. Since the total entropy is expressed 
as 



S = J Sf^J—gs'y = ^3 / 



the normalization of W should be taken to be 



for the parametrization Eq.(Bs) and 



J d 3 rW(r-f") = 1, 



-gW(f-f') = 1, 



(45) 



(46) 



(47) 



for the parametrization Eq.(|44|). In the usual SPH calculations, it is not desirable to introduce in W the space- 
time dependence through its normalization condition. In this aspect, the most natural way to introduce the SPH 
representation is Eq.(E3|). With this choice, the SPH action is given by 



I SPH = - J dr J d 3 xJ2 > 

i 



-gs~/ 



(48) 



The variational principle leads to the following equation of motion, 



dr 



1 Pi + Qi 



3_ Pj + Qj 



vt Pt + Qi ( 1 



Vi ( P + Q + e 



i 

(V500 - ufVgVi). , 



(49) 
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where 



'P + Q + e" 
TTj = ItVi ( ; | gv t 



(50) 



and the operator V is just the simple derivative operator with respect to the coordinate variable in use. 

F. Hyperbolic Coordinates 

For ultrarelativistic heavy-ion collisions, a useful set of variables are 



= V^ 2 



1 ,t + z 

ri = - tanh , 

' 2 t-z 



r T 



(51) 
(52) 

(53) 



As mentioned above, the initial conditions for RHIC processes are specified in terms of proper time rather than of 
coordinate time t. The variable r is not exactly the physical proper time of the matter, but for the initial times it 
may approximate the proper time. 

The metric tensor for this coordinate system is given by 

.9oo = 1, 

/I \ 
g= 1 , 
V r 2 / 

Since the metric is space independent, we can use the parametrization, 

n 

TUSi = s* = ^2vjW(qij) , 



3=1 



where 



and W is normalized as 



The SPH equation becomes 



4tt / q 2 dq W (q) = 1. 



d _ 1 
dr t t— 1 



1 P l + Q l 1 Pj + Q 3 



i 2 s 2 



I 2 s 2 



where the r\ component of the momentum is related to the velocity drj/dr as 



2 j P + Q + e \ dr) 



whereas in the transverse directions, we have 



7r T = wy 



P + Q + e\ dr T 
~d7 



The Lorentz factor is given by 



1 — Vj, — T 2 V 2 
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III. EXAMPLES 



We have formulated the relativistic SPH method, appropriate to the study of RHIC processes. In order to check 
its validity and efficiency we apply, in this section, our method to several known problems and compare the results to 
the analytic or known numerical solutions. 



A. Adiabatic Flow of Massless Pion Gas 



For adiabatic flows of perfect gas, the entropy is conserved. Thus Vi = and we have Q = 0. In the following 
examples, we consider the cases where Q — 0. 



1. Landau Model 



Let us first investigate a well-known analytically soluble flow. The Landau model is one of the few examples of 
relativistic fluid flow for which analytical solution is available. Consider a one-dimensional, relativistic, massless, 
baryon-free gas initially at rest. The equation of state is 

P = -e = Cs A /\ 
3 



where 

i :i 

c 



\128n 2 J 



Since we are dealing with a perfect gas, Q = 0. To apply the SPH method, we introduce the discrete one dimensional 
space variable {r/i (t) , i = 1, .., n}. The relation between the momentum and velocity is then 

7r = CW^V 73 ^ (54) 

In this case, v can be solved analytically in terms of 7r. 

In Fig.l-a and -b, we show the results of our SPH calculations, together with the exact solution M. 

Fig. 1-a, Fig.l-b 

In this example, we took only 100 particles with equally spaced {rji}. As we see, in spite of a rather small number of 
particles, the SPH solution is quite satisfactory in this example. In particular, when we use the r\ — r coordinates with 
an appropriate distribution of fj's (Fig. 1-b), an excellent agreement with the analytical solution can be obtained. 



2. 3D Scaling Solution 

A simple analytical solution for a 3 dimensional relativistic pion gas is available. It is just a generalization of the 
one-dimensional scaling solution and given by 

so 

s = 

V t 2 — % 2 — y 2 

To see the efficiency of the SPH approach presented here, we reproduce this solution numerically in the full 3D 
numerical code, without making use of the spherical symmetry. In Fig. 2, we show the result of such calculations. 

Fig.2 

As we see, our numerical calculations reproduces the analytical solution fairly well. One of the advantages of the 
SPH approach is that the coding for 3D cases is almost the same as for the ID case. The only problem is a rapid 
increases of the number of particles for higher dimensions if we want to keep a high accuracy. A direct coding would 
require a computational time proportional to n 2 , where n is the number of particles. However, in the absence of 
long range forces such as the gravitational or Coulomb interactions, we may apply techniques such as the linked-list 
method to reduce the computational time to the order of nlogn. In the example shown above, we used 50 x 50 x 50 
particles, which took less than 2 minutes for one time step in a reasonable PC(Pcntium-II). Here, we put a rather 
large number of particles, since it was a somewhat stringent test due to the divergent nature of the solution at the 
border. 
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3. Transverse Expansion on Longitudinal Scaling Expansion 

As a further test, closer to a realistic situation than that of Fig. 2, we calculated the transverse expansion of a 
cylindrically symmetric homogeneous massless pion gas, undergoing a longitudinal scaling expansion, and initially 
at rest in transverse directions. In Fig. 3, we show an example. Such a problem has been discussed by several 
authors as a useful base to understand the transverse expansion. Here, we compare our results (again a full 3D 
calculation without assuming cylindrical symmetry) with (2+1) numerical results, obtained by the use of the method 
of characteristics. Jnj] 

Fig. 3 

In this example, we used also 50 x 50 x 50 particles. The result is quite satisfactory. If we may decrease the accuracy 
by 10%, we can reduce the particle number almost by one order of magnitude. 



B. Non Adiabatic Case: Shock Waves and Artificial Viscosity 



As seen in the previous examples, our entropy-based relativistic SPH method works quite well for the adiabatic 
dynamics (Q = 0) of the massless pion gas. However, for the application to realistic problems, it is fundamental to 
see how this scheme works for non-adiabatic cases (Q ^ 0), too. For this purpose, we study some examples of one 
dimensional shock problems. 



1. Compression Shock (Normal Matter) 

Whenever there exists a shock wave, always exists the production of entropy through the shock front. The shock front 
manifests as a discontinuity in thermodynamical quantities in a hydrodynamic solution. Mathematically speaking, 
the shock front should be treated as a boundary to connect two distinct hydrodynamic solutions. To reproduce such 
a discontinuous behavior, the full degrees of freedom of hydrodynamics are required. The smoothed particle ansatz 
excludes such a possibility from the beginning. Since there do not exist short-wavelength excitation modes in the 
SPH ansatz, the energy and momentum conservation required by the hydrodynamics results in very rapidly oscillating 
motion of each particle. Such a situation occurs, for example, when a very high energy density gas is released into a 
low density region. This kind of shock, for the case of a baryon gas, is discussed in Ref.( Hj) and also, in the SPH 
context, in Ref.( fj~|]). Here, we applied our entropy-based SPH approach to the massless pion gas. 

Fig. 4 gives a typical behavior of SPH solution for such a situation, if entropy production is not taken into account. As 
discussed above, there appear rapid oscillations in thermodynamical quantities just behind the shock front. Actually, 
this oscillation appears always in numerical approaches if entropy production is not included. 

In order to avoid these oscillations, von Neuman and Richtmeyer [ pl| introduced the concept of pseudoviscosity. 
The idea is to set the dissipative pressure where the shock wave discontinuity is present. To do this, Neuman and 
Richtmeyer proposed the ansatz 

(aAxfp [p/pf, P>0 
0, p < 

for nonrelativistic one-dimensional hydrodynamics. Here, p is the mass density, Ax is the space grid size and a is a 
constant of the order of unity. In order to generalize the above pseudoviscosity for relativistic SPH case, we replace 
the quantity p/p by —9 = —d^u^ and Ax by h, where h is as before the width of the smoothing kernel W. More 
precisely, we take the following form which is a slightly modified expression suggested by Ref.( |12|), 



P 



-ah0 + (3(h0) 2 , 9<0, (55) 
0, 9>0. 



As mentioned before, Q is equivalent to the bulk viscosity and therefore there is no heat flow associated with it. 
What this artificial viscosity does is to convert the collective flow energy into the microscopic thermal energy. As a 
consequence, the total energy, that is, the sum of the collective flow energy and the internal thermal energy is still 
conserved. 

Fig. 5 is the solution of the same problem of Fig. 4, but with the entropy production taken into account. In this 
calculation, the parameters have been chosen as 

a = 2, (3 = 4 
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and h = 0.5 fm for 1000 SPH particles. As we see, the rapid oscillations have been smoothed out (and in turn, the 
numerical calculation became much more efficient). 

It is known that the overall energy- and momentum- flux conservation relates the ratio S2I s\ of entropy densities 
after and before the shock to the velocity v s of the shock front as (Hugoniot-Rankine relation) 

S2 2 (% s 2 -l) 1/4 

-rwr^^-^K- (56) 

In Fig. 6, we show the velocity of the shock front obtained in our SPH calculations as function of the entropy 
ratio(dots). Each point corresponds to the different initial condition. They are compared with the Hugoniot-Rankine 
relation Eq.(^6|) (curve). The accordance shows that our SPH calculation reproduces faithfully the conservation of 
kinetic energy and momentum of the flow through the shock front. 



2. Rarefaction Shock 



When the fluid presents a first order phase transition, a discontinuity appears in an expansion regime. This kind 
of shock wave has been discussed in connection to the QGP-hadron phase transition p^jl^|. In the present example, 
we use a simple MIT bag-model equation of state, 



P 




T <T C 

T > T r 



being B the bag constant and T c the critical temperature. 

Because decompression occurs with a constant pressure, we should have negative Q values for 8 > 0. We can use a 
similar expression for Q as before, 



Q 



0, 

p c e [(e - E h ) {s QG p - e)] \-aho - p (he) 



< 0, 
>o, 



where represents the Heaviside step function to let the viscosity be effective only in the transition region. In Fig. 
6, we show a result of our calculation and compare it to the analytic solution. The real sharp shock front is a little 
bit smoothed, but this is due to the small number of SPH particles (N — 1000 in this example) and hence a large h. 
Using the pseudoviscosity, the shock width turns out to be a few times h. 



IV. DISCUSSION AND PERSPECTIVES 



In the usual hydrodynamic computations using space grids, the symmetry of the problem is often a crucial factor 
to perform a calculation of reasonable size. The SPH method cures this aspect and furnishes a robust algorithm 
particularly appropriate to the description of processes where rapid expansions of the fluid should be treated. In 
this paper, we formulated an entropy-based SPH description of the relativistic hydrodynamics. We have shown that 
this approach is very promising for the study of ultrarelativistic nucleus-nucleus collision processes. The equations of 
motion are derived by a variational procedure from the SPH model action with respect to the Lagrangian comoving 
coordinates. This guarantees that the method furnishes the maximal efficiency for a given number of degrees of 
freedom, keeping strictly the energy and momentum conservation. For this reason, solutions can be obtained with a 
very reasonable precision, with a relatively small number of SPH particles. This is the basic advantage of the present 
method, when we analyze the event-by-event dynamics of the relativistic heavy-ion collisions. 

On the other hand, the precision of this method increases rather slowly with the number of SPH particles. Therefore, 
a relatively large number of particles is required if one wants a very precise numerical solution. However, for the 
application to the RHIC physics, we may need rather crude precision especially if we consider the dubious validity of 
the rigorous hydrodynamics. For a calculation with typically 10% errors, the SPH algorithm presented here furnishes 
a very efficient tool to study the flow phenomena in the RHIC physics. 

A fundamental difficulty of the relativistic hydrodynamics for viscous fluid p2| , p3| is that the dissipation term 
causes an intrinsic instability to the system described by Eq.(^3|). This instability basically comes from the fact 
that the dissipation term contains 6 (see Eqs.( p4| , [33| )), so that it necessarily introduces the third time-derivative 
into the equation. This means that we have to specify, at least, a part of the acceleration as the initial condition. 
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Even we specify the initial acceleration, the requirement of the internal self-consistency among the equations above 
leads to intrinsically unstable solutions. Israel proposed [^,^3) to cure these difficulties by introducing higher-order 
thermodynamics with respect to deviations from the equilibrium. In the examples presented in the present paper, 
we did not address this question and simply estimated the quantity 9 from the quantities one time step before. In 
practice, this will cause no numerical instability and the behavior of the solution is quite satisfactory. 

For the future application of the present program, we need to specify more realistic initial conditions and also to 
relate the final state to the physical observable quantities, such as particle spectra. The first point is now being in 
progress by introducing the initial energy and momentum distributions of the matter using the NEXUS algorithm. 
As for the second point, that is, the problem of particle production, the possibility of incorporation of the continuous 
emission mechanism |24| is being studied. 

This work was supported in part by PRONEX (contract no. 41.96.0886.00), FAPESP (contract nos. 98/02249-4 
and 98/00317-2), FAPERJ (contract no.E-26/150.942/99) and CNPq-Brasil. 
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FIGURE CAPTIONS 

Fig. 1: a) Entropy profiles of Landau Model for different times. The numerical solutions correspond to the equal v 
for all the particles and the Cartesian coordinate system was used, b) The same but the hyperbolic coordinate 
system was used (see text). 

Fig. 2: Three dimensional scaling solution for the massless pion gas. Cartesian coordinates are used for the transverse 
direction and the hyperbolic coordinates (77 — r) for the longitudial direction. Despite the spherical symmetry, 
the SPH calculation has been carried out in the full 3D code. 

Fig. 3: Transverse temperature profiles of a cilyndrically symmetric flow with longitudinally scaling expansion. The 
SPH results (circles) are compared with the numerical solution of the space-fixed grid method. The SPH 
calculation has been done in the full 3D code. 

Fig. 4: Shock wave in one dimensional pion gas. No viscosity is used. 

Fig. 5: After the introduction of the Q term in the SPH calculation. 

Fig. 6: Test of Hugoniot-Rankine relation. Black circles are those of the SPH calculations with different initial 
conditions, and the solid line is Eq. (|56|) . 

Fig. 7: Comparison of the SPH results with the analytical solution for the rarefaction shock. Here we took B = 
A00MeV/fm 3 . 
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